#include <functional>
namespace pl = std::placeholders;
#include <vector>
#include <algorithm>
#include <map>
#include <string>
#include <boost/array.hpp>
#include <boost/multi_array.hpp>

#include <boost/numeric/odeint.hpp>

using namespace std;
using namespace boost::numeric::odeint;

#include <iostream>
#include <fstream>
#include <cstring>
#include <sstream>




//void write_lorenz( const state_type &c , const double t )
//{
//    cout << t << '\t' << c[0] << '\t' << c[1] << '\t' << c[2] << endl;
//}

//Forward declaration для того что бы можно было хранить указатель на зерно к которому 
//относится автомат, далее описание и представление классов следуют в таком порядке, что бы 
//методы класса Automata имели доступ к методам и свойствам класса Grain
class Grain;

class Automata
{
	public:
	typedef boost::array < double, 8> state_type;
	// будем хранить состояние каждой ячейки.
	state_type c;
	state_type cbuf;

	// Идентификатор автомата.
	string uid; 
	// Идентификатор соединения, занимающего автомат
	int cid;
	// Идентификатор зерна к которому относится автомат(номер зерна в массиве)
	int gid;
	// Указатель на зерно к которому относится автомат
	Grain * grain;
	// Флаг на занятость если автомат занят соединением то isempty = 0, иначе =1
	int isempty;
	// Флаг на занятость во времени, если автомат был занят на шаге по времени то поднимается этот флаг 
	// Т.о. автомат как бы зарезервирован. 
	bool reserved;
	// Количество занятых автоматов 
	static int nonempty_automatos;
	// Массив указателей на соседей
	vector < Automata * > neighbours;
    // Количество автоматов
	static const int N = 100;
	//Массив всех автоматов
    static  Automata  pull[N][N][N];
	//TODO: Завести массив резервированных автоматов
	//static vector< vector < vector < Automata > > > pull;
	// Координаты автомата
	int i, j, k;	 
	// Конструктор 
	Automata( int, int, int);
	Automata();
	// Получить всех соседей этого автомата;
	void  get_neighbours();
	void  show_neighbours(int, int, int);
	//Смотрим у какого из соединений максимальная средняя 
	//концентрация в соседних автоматах, возвращаем номер зерна
	int  get_max_con();
	// Создаем списки соседей для всех автоматов
	static void  make_neighbours_list();
	// Возврат параметра eta  - среднего арифметического isempty по всем соседним автоматам
	double get_eta();
	void show_coords();
	//Присоединяем автомат к одному из зёрен
	void go_to_grain(Grain *);
	//Эта процедура выполняется на этапе инициализации
	void go_to_grain_init(Grain *);
	bool is_empty();
	void get_all_info();
	void make_it_beasy();
	void make_it_reserved();
	//Переписать сbuf на с 
	void update();

};

class Grain
{
	public:
	// Список uid автоматов, которые относятся к этому зерну.
	std::vector< string > auto_uids;
	// Массив указателей на автоматы, которые относятся к этому зерну.
	std::vector< Automata * > auto_pointers;
	//Конструткор
	// Идентификатор соединения, из которого формируется зерно
	double cid;
	// Номер зерна в массиве, выполняет роль идентификатора
	int grain_number;

	Grain(int,int,int,int,int,int);
	Grain();

};

//vector< vector < vector < Automata > > > Automata::pull;
// Описали статический член класса	
Automata  Automata::pull[Automata::N][Automata::N][Automata::N];
int Automata::nonempty_automatos = 0;

Automata::Automata(int ni=1, int nj=1, int nk=1): c{ {0.1 ,0.5, 0.1, 0.2, 0.1, 0, 0, 0} }
{
	
	// начальные условия по концентрациям заданы в инициализаторе ^
	//c = { 0.1 ,0.5, 0.1, 0.2, 0.1, 0, 0, 0 }; 
	//На первом предыдущее состояние аналогично текущему
	cbuf = c;
	// генерируем uid (а может просто хранить указатель на объект?)
	//uid = rand() % 10000000000;
	uid  = to_string(ni) + to_string(nj) + to_string(nk);
	//Сохранили координаты
	i = ni ; j = nj; k = nk;
	// По умолчанию автомат не занят ни чем, реакции там не протекают 
	// По мере того как автомат будет чем-нибудь занят, это значение может 
	// принимать значения 5,6,7; либо сохранять -1 в случае, если автомат свободен
	cid = -1;
	// Пока автомат не занят номер зерна будет равен -1
	gid = -1;
	// флажок на всякий случай
	isempty = 1;
	//Автомат никем не зарезервирован
	reserved = false;
	// Положили объект в массив автоматов?
	this->pull[i][j][k] = *this;


}

Automata::Automata()
{
	

}


// После всего в массиве указателей neighbours будут лежать указатели на соседей 
// данного автомата. Останется просто по ним прогуляться
void Automata::get_neighbours()
{
	// wm - walk matrix - матрица для обхода соседей
	
							//	 i   j   k			
	int wm [26][3] = { 
							// {  0, 0 , 0 },
//}
      						   {  0, 0 ,+1 },
							   {  0, 0 ,-1 },

							   {  0, +1, 0 },
							   {  0, +1,+1 },
							   {  0, +1,-1 },

							   {  0, -1, 0 },
							   {  0, -1,+1 },
							   {  0, -1,-1 },

							   { +1, 0 , 0 },
      						   { +1, 0 ,+1 },
							   { +1, 0 ,-1 },

							   { +1, +1, 0 },
							   { +1, +1,+1 },
							   { +1, +1,-1 },

							   { +1, -1, 0 },
							   { +1, -1,+1 },
							   { +1, -1,-1 },

							   { -1, 0 , 0 },
							   { -1, 0 ,+1 },
							   { -1, 0 ,-1 },
							  
							   { -1, +1, 0 },
							   { -1, +1,+1 },
							   { -1, +1,-1 },

							   { -1, -1, 0 },
							   { -1, -1,+1 },
							   { -1, -1,-1 }
							};
	int pp=0;

	for (int m = 0; m < 26 ; m++)
	{
		int ii = i + wm[m][0];
		int jj = j + wm[m][1];
		int kk = k + wm[m][2];
		if (((ii >= 0) && (ii < N) )  && (( jj >= 0) && (jj < N)) && 
		   ((kk >= 0) && ( kk < N) ) ) 
		   { 
			 //Если попадаем в границы, заносим соседа в список указателей  
			 neighbours.push_back(&pull[ii][jj][kk]);   				
		   };


	}



}

// Выводим координаты соседей автомата с координатами ii jj kk
void Automata::show_neighbours(int ii, int jj, int kk)
{

	vector < Automata * > :: iterator it;
	for (it = pull[ii][jj][kk].neighbours.begin(); it!= pull[ii][jj][kk].neighbours.end(); ++it)
		cout << (**it).i <<"	" <<(**it).j << "		"<< (**it).k << endl;
	cout << "Соседи показаны" << endl;	
	//cout << (*pull[ii][jj][kk]).neighbours.size(); 
	
}


//Возвращаем номер зерна, средняя концентрация продукта в котором максимальна 
int Automata::get_max_con()
{

	//Сюда соберем концентрации компонентов
	vector <int> con_collector(50,0);
	//TODO: разобраться с 50
	// Считаем количество соседей занятых конкретным зерном (пусть их будет не более 50)
	vector <int> grains_collector(50,0);

	vector < Automata * > :: iterator it;
	//пробежимся по соседям
	for (it = neighbours.begin(); it!= neighbours.end(); ++it)
	{
		 Automata lauto = **it;
		 if (!lauto.is_empty())
		 {
			 grains_collector[lauto.gid]++;
			 con_collector[lauto.gid] += lauto.cbuf[(**it).cid];
		 }
	}

	double max_con  = 0;
	//номер зерна с максимальной средней концентрацией
	int max_grain = -1;

	//Если хоть где-то в массиве есть один не нулевой счетчик зёрен, 
	//считаем среднии концентрации и тд, иначе всё остается по-прежнему. 
	for (int ni = 0; ni < 50; ni++ ) 
	{
		if (grains_collector[ni] != 0)
		{
			double lcon = con_collector[ni] / grains_collector[ni];
			if (max_con >= lcon) {max_con = lcon; max_grain = ni; }
		}
	}

	//возвращаем номер зерна
	return max_grain;

}

// Выводим координаты  автоматов 
void Automata::show_coords()
{

	for (int ii = 0; ii < N; ii++ ) 
		for (int jj = 0; jj < N; jj++ ) 
			for (int kk = 0; kk < N; kk++ ) 
				 cout << pull[ii][jj][kk].i << endl;

	
}



//Создаем список соседей для каждого из автоматов
void Automata::make_neighbours_list()
{
	//Формируем списки соседей
	for (int ii = 0; ii < N; ii++ ) 
		for (int jj = 0; jj < N; jj++ ) 
			for (int kk = 0; kk < N; kk++ ) 
				 pull[ii][jj][kk].get_neighbours();


}


void Automata::go_to_grain(Grain * gp)
{
	cid = gp->cid;
	make_it_reserved();
	// Положили в указатель ссылку на зерно
	grain = gp;
	gid = gp->grain_number;

}


// Для формирования зародыша
void Automata::go_to_grain_init(Grain * gp)
{
	c[2]-=0.1;
	c[(*gp).cid]+=0.1;
	make_it_beasy();
	go_to_grain(gp);
}

// Получаем значение параметра eta (характеризует количество свободных автоматов) для текущего автомата
double Automata::get_eta()
{
		
	double eta = 0; 
	vector < Automata * > :: iterator it;
	for (it = this->neighbours.begin(); it!= this->neighbours.end(); ++it)
		eta = eta + (**it).isempty;
	eta = eta / this->neighbours.size();	
	return eta;


}

// Проверяем пустой ли автомат
bool Automata::is_empty()
{
	if (this->isempty == 0) 
	{return false;}
	else 
	{return true;}

}

// Получаем номер соединения которое превалирует в окружающих автоматах
//int Automata::get_comp_num()
//{


//}

// Вся инфа по конкретному автомату 
void Automata::get_all_info()
{
	cout	<< "[" << i << " " << j << " " << k << "]: "
			<< "соед.: " << cid << " зерно: " << gid 
			<< " занят: " << isempty << " eta: " << get_eta() << "\n";
	
}

// Обозначаем что автомат занят
void Automata::make_it_beasy()
{
	isempty = 0;
	
	nonempty_automatos++;

}

// Обозначаем что автомат зарезервирован
void Automata::make_it_reserved()
{	
	reserved = true; 
}

Grain::Grain(int i, int j, int k, int delta, int cmp, int grain_num)
{	
	
	auto_uids.push_back(Automata::pull[i][j][k].uid);
	auto_pointers.push_back(&Automata::pull[i][j][k]);
	cid = cmp;
	grain_number = grain_num;
	Automata::pull[i][j][k].go_to_grain_init(this);
}


void Automata::update()
{
	cbuf = c;
}


Grain::Grain()
{	
	
}




class Phase
{
	public:
	// Вектор зёрен, принадлежащих данной фазе
	std::vector < Grain * > grains;
	// Название фазы
	string title;
	double dens;
	double mu;



};


class Mount
{
	public:
	//Свойства компонент храним в структуре
	struct components
	{
		string title;
		double dens;
		double mu;
	};


	// Количество автоматов в каждом из измерений, итого 
	// Рассматриваем кубик NxNxN
    int N = Automata::N;

	// Размер зародыша
	double delta;
	// TODO: убрать этот костыль с вынесением eta 
	// в члены класса Mount, eta должна передаваться в ф-ию ode
	double eta;

	// Описываем 3d массив автоматов
	//typedef boost::multi_array< Automata , 3> array_type;
    //typedef array_type::index index;
	// Всего у нас 8 компонент (продукты и реагенты)
	boost::array < Phase , 9 > cmp;
	// Определяем тип для вектора решения
	typedef boost::array< double , 8 > state_type;
	state_type cmid;
	// Номер граничной фазы
	int boundary_num = 8;

	// Ассоциативный Массив с зёрнами
	std::map < int, Grain > grains;
	//Счетчик зёрен
	int grains_counter = 0;

	// Времечко
	double t = 0;
	// Шаг по времени
	double dt = 1e-3;
	// Когда закончить
	double finish = 1.0;
	//Конструктор 
	Mount();
	// Задаем ОДУ которые будем решать
	void ode(const state_type &c , state_type &dcdt , double t);
	// Выводим результаты на экран
	void write_results( const state_type &c , const double t );
	// Размещение зародышей
	void put_nucleus(double delta, double cmp_num);
	// Забить всех соседей для каждого автомата;
	void get_all_neighbours();
	// Тут проходимся по всем автоматам и переключаем время
	void make_solution();
	// Тут решаем все что можно для одного автомата
	void solver(int, int, int);
	// Проверка всех автоматов на занятость
	bool check_all_automata();
	// Обходим все автоматы с различными целями, 
	// Например для того что бы обозначить занятость
	void touch_all();
	//Решатель из odeint
	runge_kutta4< state_type > rk;
	void get_results();


};


Mount::Mount()
{
	// Создаем 3Д массив автоматов
   // array_type automats(boost::extents[N][N][N]);
	
	//Инициаллизируем автоматы
	for (int i = 0; i < N; i++ ) 
		for (int j = 0; j < N; j++ ) 
			for (int k = 0; k < N; k++ ) 
				 Automata(i,j,k);
	cout << "Массив автоматов создан" << endl;
	cout << "Формируем списки соседей" << endl;
	
	Automata::pull[1][1][1].make_neighbours_list();
	
	//Формируем списки соседей
	//for (int i = 0; i < N; i++ ) 
	//	for (int j = 0; j < N; j++ ) 
	//		for (int k = 0; k < N; k++ ) 
	//			 Automata::pull[i][j][k].get_neighbours();
	
	cout << "Списки соседей сформированы" << endl;

	//(*Automata::pull[0][0][0]).show_neighbours(0,0,0);



	// Инициализация компонентов системы
	cmp[0].title = "Ti";		cmp[0].dens = 4.5400;		cmp[0].mu = 47.88;
	cmp[1].title = "N";			cmp[1].dens = 0.0012;		cmp[1].mu = 14.00;
	cmp[2].title = "Al";		cmp[2].dens = 2.6989;		cmp[2].mu = 26.98;
	cmp[3].title = "O"; 		cmp[3].dens = 0.0014;		cmp[3].mu = 16.00;
	cmp[4].title = "Si";		cmp[4].dens = 2.3300; 		cmp[4].mu = 28.08;
	cmp[5].title = "TiN"; 		cmp[5].dens = 5.4400; 		cmp[5].mu = 61.88;
	cmp[6].title = "Al2O3";		cmp[6].dens = 3.9900;		cmp[6].mu = 101.96;
	cmp[7].title = "Si3N4"; 	cmp[7].dens = 3.4400;		cmp[7].mu = 140.24;
	cmp[boundary_num].title = "boundary";  

	// Зародыши размера 5x5x5
	delta = 5;

	//Размещаем зародыши 
	put_nucleus(5,5);	put_nucleus(5,6);	
	put_nucleus(5,7);	put_nucleus(5,6);
	put_nucleus(5,5);  	put_nucleus(5,7);
	put_nucleus(5,6);	put_nucleus(5,5);
	
	
	cout << "Зародыши размещены" << endl;

	/*
	//double eta = 1.0;

		for (int i = 0; i < N; i++ ) 
		 for (int j = 0; j < N; j++ ) 
			for (int k = 0; k < N; k++ ) 
	
				 Automata::pull[i][j][k].get_all_info(); 

	*/

	cout << "Solution started..." << endl;
	make_solution();






}


//Делаем всё что надо со всеми автоматами, обычно после очередного шага по времени
void Mount::touch_all()
{
	
	for (int ll = 0; ll < 8; ll++)
				  {
					  cmid[ll] = 0; 
				  }

	
	for (int ii = 0; ii < N; ii++ ) 
		for (int jj = 0; jj < N; jj++ ) 
			for (int kk = 0; kk < N; kk++ ) 
			{
				 
				   
				  //считаем средние концентрации 
				  for (int ll = 0; ll < 8; ll++)
				  {
					  cmid[ll] = cmid[ll] + Automata::pull[ii][jj][kk].c[ll]; 
				  }

				  //TODO: переписать сbuf
				  Automata::pull[ii][jj][kk].update(); 
				 //Если автомат зарезервирован, то теперь он занят
				 if (Automata::pull[ii][jj][kk].reserved == true and Automata::pull[ii][jj][kk].is_empty() == true) 
				 {
					 Automata::pull[ii][jj][kk].make_it_beasy();
					// cout << ii << "	" << jj << "	" << kk << "	" << Automata::pull[ii][jj][kk].cid << endl;

				 }
			}

	  double total = 0;
	  cout << "Средняя концентрация: " << endl;
	  for (int ll = 0; ll < 8; ll++)
				  {
					  cmid[ll] = cmid[ll] / (N*N*N);
					  total += cmid[ll];
					  cout << ll << ":" << cmid[ll] << endl;
				  }
	  cout <<  "Итого:" << total << endl;			  


}

//Проверяем все автоматы на занятость..если есть хотя бы один свободный, возвращаем true 
bool Mount::check_all_automata()
{
	 	
		
		if (Automata::pull[1][1][1].nonempty_automatos == N*N*N) {return false;} 
		
		cout << Automata::pull[1][1][1].nonempty_automatos << endl; 
		return true;

}

void Mount::put_nucleus(double delta, double cmp_num)
{
	    //разместили зародыш и кинули его в список зерен для всей подложки
	 
		Grain new_grain = Grain( rand() % N, rand() % N, rand() % N, delta, cmp_num, grains_counter);
		grains.insert(std::pair<int, Grain>(grains_counter, new_grain ));
		//cmp[cmp_num].grains.push_back(&grains[grains_counter]);	
		//увеличили счётчик
		grains_counter++;
}




void Mount::ode( const state_type &c , state_type &dcdt , double t)
{
    	
	double w = 1.0 / ( c[0] * cmp[0].dens +  c[1] * cmp[1].dens + c[2] * cmp[2].dens + c[3] * cmp[3].dens + c[4] * cmp[4].dens +
					   c[5] * cmp[5].dens + c[6] * cmp[6].dens + c[7] * cmp[7].dens);
	double phi1 = 300 * c[0] * c[1];
	double phi2 = 800 * c[2] * c[2] * c[3] * c[3] * c[3];
	double phi3 = 1500 * c[4] * c[4] * c[4] * c[1] * c[1] * c[1] * c[1];
	
	dcdt[0] = - w * eta * cmp[0].mu * phi1; 
	dcdt[1] = - w * eta * cmp[1].mu * (phi1 + 4 * phi3); 
	dcdt[2] = - w * eta * cmp[2].mu * 2 * phi2;
	dcdt[3] = - w * eta * cmp[3].mu * 3 * phi2;
	dcdt[4] = - w * eta * cmp[4].mu * 3 * phi3;
	dcdt[5] =  w * eta * cmp[5].mu * phi1;
	dcdt[6] =  w * eta * cmp[6].mu * phi2;
	dcdt[7] =  w * eta * cmp[7].mu * phi3;
}




void Mount::make_solution()
{
	
	while (check_all_automata())
	{
		//Обходим все-все автоматы
		for (int ii = 0; ii < N; ii++ ) 
			for (int jj = 0; jj < N; jj++ ) 
				for (int kk = 0; kk < N; kk++ ) 
					{
					  	solver(ii,jj,kk);
	
					}	
		
		t = t + dt;
		touch_all();

	}

	cout << "Total time: " << t << endl;
	cout << "Write results:" << endl;
	get_results();
}

void Mount::solver(int ii, int  jj, int  kk)
{
	//Разборки с конкретным автоматом
	//Если он пустой, смотрим, можно ли его занять 
	//Чем-нибудь, если занимаем, либо если он занят, 
	//Реашем ОДУ в этом автомате, если он не занят, 
	//и занять его нельзя, выходим.
	
	
	if ( Automata::pull[ii][jj][kk].reserved == false ) 
	//Пустой
	{
		//Типа занимаем

		//Возвращаем номер зерна 
		int gid = Automata::pull[ii][jj][kk].get_max_con();
		//Если вернулся номер зерна, 
		if ( gid >= 0 ) 
			//Присоединяем автомат к нему
			{				
				//Grain *lgrain  = &grains[gid];
				Automata::pull[ii][jj][kk].go_to_grain( &grains[gid] );

			}

	} else if ( Automata::pull[ii][jj][kk].is_empty() == false )
	{
     	

		// Перед вызовом решателя возвращаем значение параметра eta для конкретного автомата
        eta = Automata::pull[ii][jj][kk].get_eta();
		//Хак отсюда 
		//http://stackoverflow.com/questions/10976078/using-boostnumericodeint-inside-the-class#new-answer?newreg=ea580eee94e648a2a840672d931b664c
		rk.do_step(std::bind(&Mount::ode, std::ref(*this) , pl::_1 , pl::_2 , pl::_3 ) , Automata::pull[ii][jj][kk].c, t , dt);
	}



}

void Mount::get_results()
{
		
		
		for (int ll = 5; ll < 8; ll++)
		{
			string fname = cmp[ll].title + ".dat";	
			ofstream fout(fname.c_str());
	        fout.precision(8);
    	    fout.setf(ios::showpoint);
			for (int ii = 0; ii < N; ii++ ) 
		    	for (int jj = 0; jj < N; jj++ ) 
					for (int kk = 0; kk < N; kk++ ) 
						{
						  if (Automata::pull[ii][jj][kk].cid == ll)
							  fout << ii / 100.0 << "	" << jj / 100.0  << "	" << kk / 100.0 << "	" 
							       << Automata::pull[ii][jj][kk].c[Automata::pull[ii][jj][kk].cid ]  << endl;

						}


			fout.close();

		}		



}



int main(int argc, char **argv)
{

	cout << "Инициализация..." << endl;
	Mount System;
	

}
